# **************************************
# --*-- coding: utf-8 --*--
# @Author  : white
# @FileName: OLS回归.py
# @Time    : 2025-08-13
# **************************************
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm

"""
OLS回归，全称是普通最小二乘法（Ordinary Least Squares）回归，
是一种最常用且基础的线性回归方法，用于估计自变量与因变量之间的线性关系。
其基本思想是：在所有可能的直线（或超平面）中，
寻找一条能够使所有数据点到该直线的垂直距离（即残差，实际值与预测值之差）的平方和最小的那条线，
作为对数据的最佳拟合线。通过最小化残差平方和，OLS能够得到自变量系数的估计值，
这些估计在一定假设下（如误差项零均值、同方差、无自相关、线性关系等）具有良好的统计性质，比如无偏性、有效性以及一致性。
在这段代码中，OLS回归用于拟合一个简单线性模型 y = β₀ + β₁ * x + ε，其中β₀是截距，β₁是x的斜率，ε是随机误差。
通过statsmodels库可以很方便地构建模型、估计参数、查看统计检验结果并作图，是一种在数据分析与统计建模中非常基础和重要的工具。

调用 sm.OLS(y, X)时，实际上是在使用 Python 的 statsmodels库来拟合一个普通最小二乘（Ordinary Least Squares, OLS）线性回归模型，
目的是基于给定的自变量矩阵 X和因变量向量 y，找出能最小化残差平方和的回归系数，
也就是找到一组参数，使得模型预测值与真实观测值之间的差的平方和达到最小，从而建立变量之间的线性关系。
在背后，sm.OLS(y, X)并不是直接进行运算，而是先创建了一个OLS 模型对象，
这个对象封装了数据（y 和 X）、模型设定（比如是否包含常数项，不过在这里你已经通过 sm.add_constant(x)把常数项加进去了）、
以及后续进行统计推断和估计的方法。当你调用 .fit()方法时，
才真正开始执行回归系数的估计过程，计算出最优的参数，同时还会计算一系列统计量，
比如标准误、t 值、p 值、R²、F 统计量等，用于评估模型的拟合效果和参数的显著性。
"""


def library():
    # 设定样本数量为100
    nsample = 100
    # 生成了一个从0到10均匀分布的100个数据点作为自变量x。
    x = np.linspace(0, 10, nsample)
    # 为了进行线性回归建模，使用sm.add_constant(x)在x前面添加了一列常数1，形成设计矩阵X，这是为了在模型中包含截距项（常数项）。
    X = sm.add_constant(x)
    # 回归的真实系数被设定为beta = [1, 10]，其中第一个元素1是截距，第二个元素10是x的斜率。
    beta = np.array([1, 10])  # 实际的维度是一个2行1列的矩阵
    # 然后通过np.random.normal(size=nsample)生成了100个服从标准正态分布的随机误差项e，代表模型中未能解释的随机扰动。
    e = np.random.normal(size=nsample)
    # 真实的因变量y是通过线性组合X与beta再加上误差e计算得到的，
    # 即 y = X * beta + e，这是典型的线性回归模型形式：y = 截距 + 斜率 * x + 误差。
    y = np.dot(X, beta) + e
    # sm.OLS(y, X)创建了一个普通最小二乘回归模型对象，
    # 其中y是因变量，X是包含常数项的自变量矩阵。
    model = sm.OLS(y, X)
    # 调用model.fit()对模型进行拟合，得到拟合结果对象results
    results = model.fit()
    # 查看拟合得到的回归系数，包括截距和x的斜率，它们应接近真实值[1, 10]，但由于误差e的存在，可能不完全相等。
    print(results.params)
    # results.summary()则输出了一份详细的回归分析报告，
    # 包括各个系数的估计值、标准误、t统计量、p值、置信区间、模型的拟合优度（如R-squared）、F统计量等统计信息，用于评估模型的显著性与解释能力。
    print(results.summary())
    # 为了直观展示回归效果，代码计算了拟合的y值（即根据回归模型预测的y值）y_fitted = results.fittedvalues，
    y_fitted = results.fittedvalues
    fig, ax = plt.subplots(figsize=(8, 6))
    # 原始数据点用圆圈'o'表示，拟合的回归线用红色虚线点线'r--.'表示，
    ax.plot(x, y, 'o', label='data')
    ax.plot(x, y_fitted, 'r--.', label='OLS')
    ax.legend(loc='best')
    # ax.axis((-0.05, 2, -1, 25))
    plt.show()


def main():
    nsample = 100
    x = np.linspace(0, 10, nsample)
    X = sm.add_constant(x)
    beta = np.array([1, 10])
    e = np.random.normal(size=nsample)
    y = np.dot(X, beta) + e
    beta_hat = np.linalg.solve(X.T @ X, X.T @ y)
    print(beta_hat) # [ 1.10839559 10.0132127 ]


if __name__ == '__main__':
    # main()
    library()
